######For the chlorophyll plot
setwd('/media/Iomega_HDD/boulot/Article/kd/last_23_01_2012/data')
load("t3.Rdata")
t3 <- as.data.frame(t3)
titick<-c(seq(0.01,0.1, by=0.01),seq(0.1,1, by=0.1),seq(1,10, by=1),seq(10,100, by=10))

opar<-par()
postscript('../Figures/Figure2.ps',width=10, height=20, horizontal=F)
# png('../Figures/Figure2.png')
#png('~/Desktop/kd_model/Paper/Figures/Figure2.png',width=1200, height=800)
par(mar=c(5,6,0.5,0.5))
par(mfrow=c(2,1))
plot(t3$gsm_chl ~ t3$chla_avg, log="xy", xlim=c(0.01,100), ylim=c(0.01,100), pch=16, cex=1.2,
 xlab=expression(paste(Chl[in_situ], "  (mg ", m^-3,")")), ylab=expression(paste(Chl[GSM], "  (mg ", m^-3,")")), cex.axis=2, cex.lab=2, xaxt="n", yaxt="n")
axis(side=1, at=titick, labels=F)
axis(side=2, at=titick, labels=F)
axis(side=1, at=c(0.01,0.1,1,10,100), labels=c("0.01","0.1","1","10","100"), cex.axis=2)
axis(side=2, at=c(0.01,0.1,1,10,100), , labels=c("0.01","0.1","1","10","100"), cex.axis=2)
abline(a=0,b=1,lwd=2, col="dark gray" )
model<-lm(log(t3$gsm_chl) ~ log(t3$chla_avg))
#abline(model, lwd=2, col="dark red")
curve(0.5*x,0.0001,200, col="orange", lwd=3, lty=2, add=T)
curve(1.5*x,0.0001,200, col="orange", lwd=3, lty=2, add=T)
text(0.01,75,pos=4,expression(paste(R^2," = ", sep="")), cex=2)
text(0.03,75,pos=4,paste(round(summary(model)$adj.r.square,3), sep=""), cex=2)
text(0.01,35,pos=4,paste("slope = ",round(summary(model)$coefficients[2,1],3)), cex=2)
text(0.01,18,pos=4,paste("N= ",length(which((t3$gsm_chl)> -100)),sep=""), cex=2)
text(100, 0.10, "(a)", cex=2)
legend("bottomright", lty=c(1,2), lwd=c(2,2), col=c("dark gray", "orange"), legend=c("1-to-1 line","50 percent confidence interval"), cex=1.5, bty="n")


######For the Acdm plot
titick2<-c(seq(0.001,0.01, by=0.001),seq(0.01,0.1, by=0.01),seq(0.1,1, by=0.1),seq(1,10, by=1))
cdm443 <- t3$ag443+t3$ad443
     
plot(t3$gsm_cdm~ cdm443, log="xy", xlim=c(0.001,10), ylim=c(0.001,10), pch=16, cex=1.2,
 xlab=expression(paste(a[g](443)," + ",a[d](443), "  (", m^-1,")")), ylab=expression(paste(a[CDM]^GSM, "  (", m^-1,")")), cex.axis=2, cex.lab=2, xaxt="n", yaxt="n")
abline(a=0,b=1,lwd=3, col="dark gray" )
axis(side=1, at=titick2, labels=F)
axis(side=2, at=titick2, labels=F)
axis(side=1, at=c(0.001,0.01,0.1,1,10), labels=c("0.001","0.01","0.1","1","10"), cex.axis=2)
axis(side=2, at=c(0.001,0.01,0.1,1,10), , labels=c("0.001","0.01","0.1","1","10"), cex.axis=2)
model<-lm(log(t3$gsm_cdm)~ log(cdm443))
#abline(model, lwd=2, col="dark red")
curve(0.5*x,0.0001,200, col="orange", lwd=3, lty=2, add=T)
curve(1.5*x,0.0001,200, col="orange", lwd=3, lty=2, add=T)
text(0.001,6.5,pos=4,expression(paste(R^2," = ", sep=" ")), cex=2)
text(0.004,6.5, pos=4, paste(round(summary(model)$adj.r.square,3)), cex=2)     
text(0.001,3.25,pos=4,paste("slope = ",round(summary(model)$coefficients[2,1],3), sep=""), cex=2)
text(0.001,1.7,pos=4,paste("N = ",length(which((cdm443)> -100)), sep=""), cex=2)
text(10, 0.010, "(b)", cex=2)
legend("bottomright", lty=c(1,2), lwd=c(2,2), col=c("dark gray", "orange"), legend=c("1-to-1 line","50 percent confidence interval"), cex=1.5, bty="n")

dev.off()
par(opar)


######For the bbp plot
titick2<-c(seq(0.001,0.01, by=0.001),seq(0.01,0.1, by=0.01),seq(0.1,1, by=0.1),seq(1,10, by=1))

opar<-par()
postscript('../Figures/Figure2bis.ps',width=10, height=20, horizontal=F)
# png('../Figures/Figure2bis.png')
par(mar=c(5,6,0.5,0.5))
plot(t3$gsm_bbp[which(t3$bb_qc==0)]~ t3$bb443[which(t3$bb_qc==0)], log="xy", xlim=c(0.0001,10), ylim=c(0.0001,10), pch=16, cex=1.2,
xlab=expression(paste(b[b443]^"in situ", sep="")), ylab=expression(paste(b[b(443)]^GSM)), cex.axis=2, 
     cex.lab=2, xaxt="n", yaxt="n")
abline(a=0,b=1,lwd=3, col="dark gray" )
axis(side=1, at=titick2, labels=F)
axis(side=2, at=titick2, labels=F)
axis(side=1, at=c(0.001,0.01,0.1,1,10), labels=c("0.001","0.01","0.1","1","10"), cex.axis=2)
axis(side=2, at=c(0.001,0.01,0.1,1,10), , labels=c("0.001","0.01","0.1","1","10"), cex.axis=2)
model<-lm(log(t3$gsm_bbp[which(t3$bb_qc==0)])~ log(t3$bb443[which(t3$bb_qc==0)]))
#abline(model, lwd=2, col="dark red")
curve(0.5*x,0.0001,200, col="orange", lwd=3, lty=2, add=T)
curve(1.5*x,0.0001,200, col="orange", lwd=3, lty=2, add=T)
text(0.0001,4,pos=4,expression(paste(R^2,"= ", sep=" ")), cex=2)
text(0.0003,4, pos=4, paste(round(summary(model)$adj.r.square,3)), cex=2)     
text(0.0001,0.8,pos=4,paste("slope = ",round(summary(model)$coefficients[2,1],3), sep=""), cex=2)
text(0.0001,0.1,pos=4,paste("N = ",length(which((t3$bb443)> -100)), sep=""), cex=2)
legend("bottomright", lty=c(1,2), lwd=c(2,2), col=c("dark gray", "orange"), legend=c("1-to-1 line","50 percent confidence interval"), cex=1.5, bty="n")

dev.off()
par(opar)

